function handles = compareGroups(handles)
% compareGroups
%
% Editing:
% gs may 2005
% ----------------------------------------------------------

% figure
hfig = figure;
set(hfig, 'Units', 'normalized', 'Position', [.15 .15 .7 .7], 'DoubleBuffer', 'on', 'Color', 'w');

subCols = 4;
subRows = 3;
plotnum = 0;

% =========================== AVERAGE MAPS ===============================================

for Ngroup = 1:2
    
    % process group A then group B
    if Ngroup == 1
        set(handles.groupList, 'Value', handles.analysis.mapAverager.groupA.groupVal);
    elseif Ngroup == 2
        set(handles.groupList, 'Value', handles.analysis.mapAverager.groupB.groupVal);
    end
    handles = updateLists(handles);
    handles = updateCurrentGroup(handles);
    
    % get list of names of cells/expts
    listOfCells = get(handles.cellList, 'String');
    
    % alignment options
    alignOptions = 'alignBySliceXY';
    %     alignOptions = 'alignBySomaXY';
    %     alignOptions = 'alignBySomaXSliceY';
    % TODO: incorporate these options into the gui
    
    % main loading loop
    numCells = size(listOfCells, 1);
    mapSpacing = [];
    cellType = [];
    XY = [];
    barrelXY = [];
    X = [];
    Y = [];
    M = [];
    
    % type of map to analyze
    handles = getMapStr(handles);
    mapStr = handles.analysis.mapAverager.mapType;

    for n = 1:numCells
        
        % load the cells' data
        cellA = eval([listOfCells{n}]);     % load the cell/expt data
        mapSpacing(n) = cellA.mapSpacing;   % map spacing
        try
            cellType{n} = cellA.cellType;       % cell type
        catch
            cellType{n} = 'neuron';
        end
        if isfield(cellA, 'xyNew')
            XY(:,:,n) = [cellA.xyNew(1) cellA.xyNew(2) * (-1)];           % soma XY position
        else
            XY(:,:,n) = cellA.xy;               % soma XY position
            disp(cellA.experiment)
        end
        try
            barrelXY(:,:,n) = cellA.barrelXY;   % barrel XY position
        catch
        end
        
        eval(['map = ' mapStr]);
%         map = cellA.genMean;                % cell's genMean map
        M(:,:,n) = map;
        
        % re-orient if necessary
        try
            if strcmp(cellA.orientation, 'upside down') % Y values need to change sign
                %         XY(:,2,n) = -XY(:,2,n);
                %         barrelXY(:,2,n) = -barrelXY(:,2,n);
                M(:,:,n) = flipud(map);
            end
        catch
        end
    end
    
    % some tidying up etc
    maps = M;
    maps(isnan(maps)) = inf;
    XY = squeeze(XY)';
    try
        barrelXY = squeeze(barrelXY)';
    catch
    end
    % X = XY(:,1) + barrelXY(:,1);
    X = XY(:,1);
    Y = XY(:,2);
    
    % distance scale
    if numel(unique(mapSpacing)) ~= 1
        beep;
        disp('Error in mapAvgFromCellList: map spacing apparently varies across cells.')
        return
    else
        [r,c] = size(maps(:,:,1));
        dx = unique(mapSpacing);
        xdata = [0:dx:(r-1)*dx] - (r-1)*dx/2;
        ydata = xdata;
    end
    
%     % shift maps to align to barrel/septal grid
%     try
%         switch alignOptions
%             case 'none'
%             case 'alignBySliceXY'
%                 for n = 1:numCells
%                     newmap(:,:,n) = cellShifter(maps(:,:,n), barrelXY(n,:), dx);
%                     X(n) = X(n) - barrelXY(n,1);
%                     Y(n) = Y(n) - barrelXY(n,2);
%                 end
%                 maps = newmap((1+r/2):(r+r/2),(1+c/2):(c+c/2),:);
%                 maps(isnan(maps)) = inf;
%             case 'alignBySomaXY'
%                 for n = 1:numCells
%                     newmap(:,:,n) = cellShifter(maps(:,:,n), XY(n,:), dx);
%                     X(n) = X(n) - XY(n,1);
%                     Y(n) = Y(n) - XY(n,2);
%                 end
%                 maps = newmap((1+r/2):(r+r/2),(1+c/2):(c+c/2),:);
%                 maps(isnan(maps)) = inf;
%             case 'alignBySomaXSliceY'
%                 for n = 1:numCells
%                     newmap(:,:,n) = cellShifter(maps(:,:,n), [XY(n,1) barrelXY(n,2)], dx);
%                     X(n) = X(n) - XY(n,1);
%                     Y(n) = Y(n) - barrelXY(n,2);
%                 end
%                 maps = newmap((1+r/2):(r+r/2),(1+c/2):(c+c/2),:);
%                 maps(isnan(maps)) = inf;
%         end
%     catch
%         error(lasterr)
%     end
%     numCells, size(maps)

    % average map
    maps(maps==inf) = NaN;
    avgMap = mapAverage(maps, 1);
    avgMap(isnan(avgMap)) = inf;
    maps(isnan(maps)) = inf;
    
    % output:
    if Ngroup == 1
        handles.analysis.mapAverager.groupA.avgMap = avgMap;
        handles.analysis.mapAverager.groupA.mapStack = maps;
        handles.analysis.mapAverager.groupA.xdata = xdata;
        handles.analysis.mapAverager.groupA.ydata = ydata;
    elseif Ngroup == 2
        handles.analysis.mapAverager.groupB.avgMap = avgMap;
        handles.analysis.mapAverager.groupB.mapStack = maps;
        handles.analysis.mapAverager.groupB.xdata = xdata;
        handles.analysis.mapAverager.groupB.ydata = ydata;
    end
    
    % plot of average map
    plotnum = plotnum + 1;
    hax(plotnum) = subplot(subRows, subCols, plotnum);
    himg = imagesc(avgMap);
    set(himg, 'XData', xdata, 'YData', ydata);
    set(hax(plotnum), 'XLim', [min(xdata)-dx/2 max(xdata)+dx/2], 'YLim', [min(ydata)-dx/2 max(ydata)+dx/2])
    pbaspect([1 1 1]);
    colormap(flipud(colormap(jet2)));
    title(['Average map (N=' num2str(numCells) ')']);
    
    % store map data
    set(himg, 'UserData', maps);
    
    % pixel/interpolate option via uicontextmenu
    cmenu = uicontextmenu;
    set(himg, 'UIContextMenu', cmenu);
    % cb1 = ['pix2interp'];
    item1 = uimenu(cmenu, 'Label', 'activate pix/interp options', 'Callback', 'pix2interp');
    pix2interp(himg);
    interpRedOutForAvg(himg, 1);
    
    % plot the soma positions 
    hold on;
    faceColor = ['b'; 'r'];
    for n = 1:numel(X)
        if strcmp(cellType{n}, 'pyramidal') | strcmp(cellType{n}, 'pyramid')
            marker = '^';
        elseif strcmp(cellType{n}, 'stellate')
            marker = 'pentagram';
        elseif strcmp(cellType{n}, 'interneuron')
            marker = 'o';
        else
            marker = 'square';
        end
        hplt = [];
        hplt(n) = plot(X(n), Y(n), 'k^');
        set(hplt(n), 'Marker', marker, 'MarkerEdgeColor', 'k', ...
            'MarkerFaceColor', faceColor(Ngroup), 'MarkerSize', 6);
    end
    somaXmean(Ngroup) = mean(X);
    somaYmean(Ngroup) = mean(Y);
    
    pos = get(gca, 'Position');
    set(gca, 'Position', [pos(1) pos(2) pos(3) pos(4)*.9]);
    set(gca, 'Tag', ['avgMap_' num2str(Ngroup)]);
    
    titleString = handles.analysis.mapAverager.currentGroup;
    ht = text(0, 0, titleString);
    set(ht, 'Units', 'normalized', 'Position', [-.2 1.2 0], ...
        'FontSize', 11, 'FontWeight', 'bold', ...
        'Interpreter', 'none');
    
%     if Ngroup == 1
%         ht = text(0, 0, 'R-click on avg maps or diff image for interpolation options');
%         set(ht, 'Units', 'normalized', 'Position', [-.6 -.6 0]);
%     end
    
%     colorbar;
    hold off;
    
    % plot of individual maps
    mapN = 1;
    plotnum = plotnum + 1;
    hax(plotnum) = subplot(subRows, subCols, plotnum);
    himg = imagesc(maps(:,:,mapN));
    set(himg, 'XData', xdata, 'YData', ydata);
    set(hax(plotnum), 'XLim', [min(xdata)-dx/2 max(xdata)+dx/2], 'YLim', [min(ydata)-dx/2 max(ydata)+dx/2])
    pbaspect([1 1 1]);
    
    % store map and soma data
    data.maps = maps;
    data.X = X;
    data.Y = Y;
    data.cellType = cellType;
    data.cellNames = listOfCells;
    % data.cellNames = filename;
    data.currentNum = 1;
    set(himg, 'UserData', data);
    
    % page through maps option via ButtonDownFcn
    set(himg, 'ButtonDownFcn', 'pageThroughMaps(gco)');
    
    % title(data.cellNames(data.currentNum, :));
    title(data.cellNames{data.currentNum});
    
    % plot the soma position
    hold on;
    n = mapN;
    if strcmp(cellType{n}, 'pyramidal') | strcmp(cellType{n}, 'pyramid')
        marker = '^';
    elseif strcmp(cellType{n}, 'stellate')
        marker = 'pentagram';
    elseif strcmp(cellType{n}, 'interneuron')
        marker = 'o';
    else
        marker = 'square';
    end
    hplt = [];
    hplt(n) = plot(X(n), Y(n), 'k^');
    set(hplt(n), 'Marker', marker, 'MarkerEdgeColor', 'k', ...
        'MarkerFaceColor', 'w', 'MarkerSize', 6, 'Tag', 'soma marker');
    
    % OPTION to turn on pixvals for measuring distances
    % pixvals on
    % pixvals off
    pos = get(gca, 'Position');
    set(gca, 'Position', [pos(1) pos(2) pos(3) pos(4)*.9]);
    
    colorbar;
    hold off;
    
end

% =========================== CLIMs ==========================================

lims = [get(hax(1), 'CLim') get(hax(3), 'CLim')];
minLim = min(lims);
maxLim = max(lims);
set(hax([1 3]), 'CLim', [minLim maxLim]);
axes(hax(1));
hcb(1) = colorbar;
axes(hax(3));
hcb(2) = colorbar;
set(hcb, 'Tag', 'avgMapCB');

cmenu = uicontextmenu;
set(hax([1 3]), 'UIContextMenu', cmenu);
item1 = uimenu(cmenu, 'Label', 'Change CLim', 'Callback', 'changeCLim');


% =========================== MAP SUBTRACTION ===============================================

% compute difference image
diffMap = handles.analysis.mapAverager.groupB.avgMap - handles.analysis.mapAverager.groupA.avgMap;

% plot of average map
plotnum = plotnum + 1;
plotnum = plotnum + 1;
hax(plotnum) = subplot(subRows, subCols, plotnum);
himg = imagesc(diffMap);
set(himg, 'XData', xdata, 'YData', ydata);
set(hax(plotnum), 'XLim', [min(xdata)-dx/2 max(xdata)+dx/2], 'YLim', [min(ydata)-dx/2 max(ydata)+dx/2])
pbaspect([1 1 1]);
colormap(flipud(colormap(jet2)));
% title(['Average map (N=' num2str(numCells) ')'])

% store map data
set(himg, 'UserData', maps);
% pixel/interpolate option via uicontextmenu
cmenu = uicontextmenu;
set(himg, 'UIContextMenu', cmenu);
% cb1 = ['pix2interp'];
item1 = uimenu(cmenu, 'Label', 'activate pix/interp options', 'Callback', 'pix2interp');
pix2interp(himg);
% red-out 
M = get(himg, 'CData');
M(M==-inf) = inf;
set(gco, 'CData', M);

% TODO:
hold on;
hplt(n) = plot(somaXmean(1), somaYmean(1), 'k^');
set(hplt(n), 'Marker', '^', 'MarkerEdgeColor', 'k', 'MarkerFaceColor', 'b', 'MarkerSize', 8);
hplt(n) = plot(somaXmean(2), somaYmean(2), 'kv');
set(hplt(n), 'Marker', '^', 'MarkerEdgeColor', 'k', 'MarkerFaceColor', 'r', 'MarkerSize', 8);

pos = get(gca, 'Position');
set(gca, 'Position', [pos(1) pos(2) pos(3) pos(4)*.9]);

titleString = ['Difference image: <' handles.analysis.mapAverager.groupB.groupName ...
        '> minus <' handles.analysis.mapAverager.groupA.groupName '>'];
ht = text(0, 0, titleString);
set(ht, 'Units', 'normalized', 'Position', [-.4 1.1 0], ...
    'FontSize', 11, 'FontWeight', 'bold', ...
    'Interpreter', 'none');

% ht = text(0, 0, 'R-click on map for interpolation options');
% set(ht, 'Units', 'normalized', 'Position', [-.2 -.1 0]);

colorbar;
hold off;

% =========================== ROI ANALYSES ===============================================

% plot of ROIs
plotnum = 9;
hax(plotnum) = subplot(subRows, subCols, plotnum);
mapEdge = size(handles.analysis.mapAverager.groupA.avgMap,1);
[X,Y] = meshgrid((1:mapEdge));
hplt = plot(X,Y);
set(hplt, 'Color', [.7 .7 .7], 'LineStyle', 'none', 'Marker', '.');
set(gca, 'Visible', 'off');
pos = get(gca, 'Position');
set(gca, 'Position', [pos(1) pos(2) pos(3) pos(4)*.9]);
hold on;
% plot col avg ROI:
rowStart = handles.analysis.mapAverager.ROI.colROIfirstRow;
rowEnd = handles.analysis.mapAverager.ROI.colROIlastRow;
colStart = handles.analysis.mapAverager.ROI.colROIfirstCol;
colEnd = handles.analysis.mapAverager.ROI.colROIlastCol;
pos = [colStart-.5 rowStart-.5 (colEnd - colStart)+1 (rowEnd - rowStart)+1];
hrect = rectangle('Position', pos);
set(hrect, 'EdgeColor', 'k', 'LineWidth', 2);
% plot row avg ROI:
rowStart = handles.analysis.mapAverager.ROI.rowROIfirstRow;
rowEnd = handles.analysis.mapAverager.ROI.rowROIlastRow;
colStart = handles.analysis.mapAverager.ROI.rowROIfirstCol;
colEnd = handles.analysis.mapAverager.ROI.rowROIlastCol;
pos = [colStart-.5 rowStart-.5 (colEnd - colStart)+1 (rowEnd - rowStart)+1];
hrect = rectangle('Position', pos);
set(hrect, 'EdgeColor', [.35 .55 0], 'LineWidth', 2);
% plot region avg ROI:
rowStart = handles.analysis.mapAverager.ROI.allROIfirstRow;
rowEnd = handles.analysis.mapAverager.ROI.allROIlastRow;
colStart = handles.analysis.mapAverager.ROI.allROIfirstCol;
colEnd = handles.analysis.mapAverager.ROI.allROIlastCol;
pos = [colStart-.5 rowStart-.5 (colEnd - colStart)+1 (rowEnd - rowStart)+1];
hrect = rectangle('Position', pos);
set(hrect, 'EdgeColor', [0 .35 .55], 'LineWidth', 2);
% tidy up
set(gca, 'YDir', 'reverse');
axis tight;
xlims = get(gca, 'XLim');
ylims = get(gca, 'YLim');
set(gca, 'XLim', [xlims(1)-1 xlims(2)+2]);
set(gca, 'YLim', [ylims(1)-1 ylims(2)+2]);
daspect([1 1 1]);

% ROI averages -- group A
set(handles.groupList, 'Value', handles.analysis.mapAverager.groupA.groupVal);
handles = updateLists(handles);
handles = updateCurrentGroup(handles);
handles = mapAvgFromCellListB(handles, 0);
% col avg
handles.analysis.mapAverager.ROI.option = 'col';
handles = analyzeROI(handles);
handles.analysis.mapAverager.groupA.colROI.vecX = handles.analysis.mapAverager.ROI.vecX;
handles.analysis.mapAverager.groupA.colROI.vecY = handles.analysis.mapAverager.ROI.vecY;
handles.analysis.mapAverager.groupA.colROI.vecAvg = handles.analysis.mapAverager.ROI.vecAvg;
handles.analysis.mapAverager.groupA.colROI.vecErr = handles.analysis.mapAverager.ROI.vecErr;
% row avg
handles.analysis.mapAverager.ROI.option = 'row';
handles = analyzeROI(handles);
handles.analysis.mapAverager.groupA.rowROI.vecX = handles.analysis.mapAverager.ROI.vecX;
handles.analysis.mapAverager.groupA.rowROI.vecY = handles.analysis.mapAverager.ROI.vecY;
handles.analysis.mapAverager.groupA.rowROI.vecAvg = handles.analysis.mapAverager.ROI.vecAvg;
handles.analysis.mapAverager.groupA.rowROI.vecErr = handles.analysis.mapAverager.ROI.vecErr;
% region avg
handles.analysis.mapAverager.ROI.option = 'all';
handles = analyzeROI(handles);
handles.analysis.mapAverager.groupA.allROI.vecX = handles.analysis.mapAverager.ROI.vecX;
handles.analysis.mapAverager.groupA.allROI.vecY = handles.analysis.mapAverager.ROI.vecY;
handles.analysis.mapAverager.groupA.allROI.vecAvg = handles.analysis.mapAverager.ROI.vecAvg;
handles.analysis.mapAverager.groupA.allROI.vecErr = handles.analysis.mapAverager.ROI.vecErr;

% ROI averages -- group B
set(handles.groupList, 'Value', handles.analysis.mapAverager.groupB.groupVal);
handles = updateLists(handles);
handles = updateCurrentGroup(handles);
handles = mapAvgFromCellListB(handles, 0);
% col avg
handles.analysis.mapAverager.ROI.option = 'col';
handles = analyzeROI(handles);
handles.analysis.mapAverager.groupB.colROI.vecX = handles.analysis.mapAverager.ROI.vecX;
handles.analysis.mapAverager.groupB.colROI.vecY = handles.analysis.mapAverager.ROI.vecY;
handles.analysis.mapAverager.groupB.colROI.vecAvg = handles.analysis.mapAverager.ROI.vecAvg;
handles.analysis.mapAverager.groupB.colROI.vecErr = handles.analysis.mapAverager.ROI.vecErr;
% row avg
handles.analysis.mapAverager.ROI.option = 'row';
handles = analyzeROI(handles);
handles.analysis.mapAverager.groupB.rowROI.vecX = handles.analysis.mapAverager.ROI.vecX;
handles.analysis.mapAverager.groupB.rowROI.vecY = handles.analysis.mapAverager.ROI.vecY;
handles.analysis.mapAverager.groupB.rowROI.vecAvg = handles.analysis.mapAverager.ROI.vecAvg;
handles.analysis.mapAverager.groupB.rowROI.vecErr = handles.analysis.mapAverager.ROI.vecErr;
% region avg
handles.analysis.mapAverager.ROI.option = 'all';
handles = analyzeROI(handles);
handles.analysis.mapAverager.groupB.allROI.vecX = handles.analysis.mapAverager.ROI.vecX;
handles.analysis.mapAverager.groupB.allROI.vecY = handles.analysis.mapAverager.ROI.vecY;
handles.analysis.mapAverager.groupB.allROI.vecAvg = handles.analysis.mapAverager.ROI.vecAvg;
handles.analysis.mapAverager.groupB.allROI.vecErr = handles.analysis.mapAverager.ROI.vecErr;

% STATISTICAL COMPARISONS
% col avg tests
H = []; P = [];
testTypeVal = get(handles.testType, 'Value');
testTypeStr = get(handles.testType, 'String');
testType = testTypeStr{testTypeVal};
for j = 1 : size(handles.analysis.mapAverager.groupA.colROI.vecY,2)
    groupAdata = handles.analysis.mapAverager.groupA.colROI.vecY(:,j);
    groupAdata(isnan(groupAdata)) = [];
    groupBdata = handles.analysis.mapAverager.groupB.colROI.vecY(:,j);
    groupBdata(isnan(groupBdata)) = [];
    eval(['[H(j),P(j)] = ' testType '(groupAdata,' ...
        'groupBdata);']);
    if strcmp(testType, 'ranksum')
        handles.analysis.mapAverager.groupB.colROI.ttest = P;
        handles.analysis.mapAverager.groupB.colROI.pval = H;
    else
        handles.analysis.mapAverager.groupB.colROI.ttest = H;
        handles.analysis.mapAverager.groupB.colROI.pval = P;
    end
end
% row avg tests
H = []; P = [];
for j = 1 : size(handles.analysis.mapAverager.groupA.rowROI.vecY,2)
    groupAdata = handles.analysis.mapAverager.groupA.rowROI.vecY(:,j);
    groupAdata(isnan(groupAdata)) = [];
    groupBdata = handles.analysis.mapAverager.groupB.rowROI.vecY(:,j);
    groupBdata(isnan(groupBdata)) = [];
    eval(['[H(j),P(j)] = ' testType '(groupAdata,' ...
        'groupBdata);']);
    if strcmp(testType, 'ranksum')
        handles.analysis.mapAverager.groupB.rowROI.ttest = P;
        handles.analysis.mapAverager.groupB.rowROI.pval = H;
    else
        handles.analysis.mapAverager.groupB.rowROI.ttest = H;
        handles.analysis.mapAverager.groupB.rowROI.pval = P;
    end
end
% region avg tests
H = []; P = [];
groupAdata = handles.analysis.mapAverager.groupA.allROI.vecY;
groupAdata(isnan(groupAdata)) = [];
groupBdata = handles.analysis.mapAverager.groupB.allROI.vecY;
groupBdata(isnan(groupBdata)) = [];
eval(['[H,P] = ' testType '(groupAdata, groupBdata);']);
if strcmp(testType, 'ranksum')
    handles.analysis.mapAverager.groupB.allROI.ttest = P;
    handles.analysis.mapAverager.groupB.allROI.pval = H;
else
    handles.analysis.mapAverager.groupB.allROI.ttest = H;
    handles.analysis.mapAverager.groupB.allROI.pval = P;
end


if get(handles.flipSignCheck, 'Value')
    S = -1;
else
    S = 1;
end

% COLUMN average -- plot the results
plotnum = plotnum + 1;
hax(plotnum) = subplot(subRows, subCols, plotnum);
pos = get(gca, 'Position');
set(gca, 'Position', [pos(1) pos(2) pos(3) pos(4)*.85]);
h = errorbar(handles.analysis.mapAverager.groupA.colROI.vecX, ...
    S * handles.analysis.mapAverager.groupA.colROI.vecAvg, ...
    S * handles.analysis.mapAverager.groupA.colROI.vecErr, 'bo-');
title('Column avg & SEM');
% titleString = ['Analysis method: ' titleTxt];
% ht = text(0, 0, titleString);
% set(ht, 'Units', 'normalized', 'Position', [.2 1.2 0], 'FontSize', 11);
hold on;
h = errorbar(handles.analysis.mapAverager.groupB.colROI.vecX, ...
    S * handles.analysis.mapAverager.groupB.colROI.vecAvg, ...
    S * handles.analysis.mapAverager.groupB.colROI.vecErr, 'ro-');
axis tight;
% color in the significant-difference points
hvals = handles.analysis.mapAverager.groupB.colROI.ttest;
I = find(hvals);
vecAvgSig = handles.analysis.mapAverager.groupB.colROI.vecAvg(I);
handles.analysis.mapAverager.groupB.colROI.vecAvgSig = vecAvgSig;
xvals = handles.analysis.mapAverager.groupB.colROI.vecX;
vecXSig = xvals(I);
handles.analysis.mapAverager.groupB.colROI.vecXSig = vecXSig;
h = plot(handles.analysis.mapAverager.groupB.colROI.vecXSig, ...
    S * handles.analysis.mapAverager.groupB.colROI.vecAvgSig, 'ro');
set(h, 'MarkerFaceColor', 'r');

% ROW average -- plot the results
plotnum = plotnum + 1;
hax(plotnum) = subplot(subRows, subCols, plotnum);
pos = get(gca, 'Position');
set(gca, 'Position', [pos(1) pos(2) pos(3) pos(4)*.85]);
h = errorbar(handles.analysis.mapAverager.groupA.rowROI.vecX, ...
    S * handles.analysis.mapAverager.groupA.rowROI.vecAvg, ...
    S * handles.analysis.mapAverager.groupA.rowROI.vecErr, 'bo-');
title('Row avg & SEM');
% titleString = ['Analysis method: ' titleTxt];
% ht = text(0, 0, titleString);
% set(ht, 'Units', 'normalized', 'Position', [.2 1.2 0], 'FontSize', 11);
hold on;
h = errorbar(handles.analysis.mapAverager.groupB.rowROI.vecX, ...
    S * handles.analysis.mapAverager.groupB.rowROI.vecAvg, ...
    S * handles.analysis.mapAverager.groupB.rowROI.vecErr, 'ro-');
axis tight;
% color in the significant-difference points
hvals = handles.analysis.mapAverager.groupB.rowROI.ttest;
I = find(hvals);
vecAvgSig = handles.analysis.mapAverager.groupB.rowROI.vecAvg(I);
handles.analysis.mapAverager.groupB.rowROI.vecAvgSig = vecAvgSig;
xvals = handles.analysis.mapAverager.groupB.rowROI.vecX;
vecXSig = xvals(I);
handles.analysis.mapAverager.groupB.rowROI.vecXSig = vecXSig;
h = plot(handles.analysis.mapAverager.groupB.rowROI.vecXSig, ...
    S * handles.analysis.mapAverager.groupB.rowROI.vecAvgSig, 'ro');
set(h, 'MarkerFaceColor', 'r');

% REGION average -- plot the results
pointColor = [.5 .5 .5];
plotnum = plotnum + 1;
hax(plotnum) = subplot(subRows, subCols, plotnum);
pos = get(gca, 'Position');
set(gca, 'Position', [pos(1) pos(2) pos(3) pos(4)*.85]);
% group A -- indiv
h = plot(handles.analysis.mapAverager.groupA.allROI.vecX, ...
    S * handles.analysis.mapAverager.groupA.allROI.vecY, 'ko');
set(h, 'Color', pointColor, 'MarkerSize', 4);
% group A -- avg
hold on;
h = errorbar(handles.analysis.mapAverager.groupA.allROI.vecX, ...
    S * handles.analysis.mapAverager.groupA.allROI.vecAvg, ...
    S * handles.analysis.mapAverager.groupA.allROI.vecErr, 'bo-');
title('Region avg & SEM');
% titleString = ['Analysis method: ' titleTxt];
% ht = text(0, 0, titleString);
% set(ht, 'Units', 'normalized', 'Position', [.2 1.2 0], 'FontSize', 11);
%
% group B -- indiv
h = plot(handles.analysis.mapAverager.groupB.allROI.vecX+1, ...
    S * handles.analysis.mapAverager.groupB.allROI.vecY, 'ko');
set(h, 'Color', pointColor, 'MarkerSize', 4);
% group B -- avg
hold on;
h = errorbar(handles.analysis.mapAverager.groupB.allROI.vecX+1, ...
    S * handles.analysis.mapAverager.groupB.allROI.vecAvg, ...
    S * handles.analysis.mapAverager.groupB.allROI.vecErr, 'ro-');
% axis tight;
set(gca, 'XLim', [.51 2.49], 'XTick', [1 2], 'XTickLabel', ['Group A'; 'Group B']); 
% xlabel('Group A        Group B');
% color in the significant-difference points
hvals = handles.analysis.mapAverager.groupB.allROI.ttest;
I = find(hvals);
vecAvgSig = handles.analysis.mapAverager.groupB.allROI.vecAvg(I);
handles.analysis.mapAverager.groupB.allROI.vecAvgSig = vecAvgSig;
xvals = handles.analysis.mapAverager.groupB.allROI.vecX;
vecXSig = xvals(I);
handles.analysis.mapAverager.groupB.allROI.vecXSig = vecXSig;
h = plot(handles.analysis.mapAverager.groupB.allROI.vecXSig+1, ...
    S * handles.analysis.mapAverager.groupB.allROI.vecAvgSig, 'ro');
set(h, 'MarkerFaceColor', 'r');
% tweak ylims
% ylims = get(gca, 'YLim');
% if S * handles.analysis.mapAverager.groupA.allROI.vecAvg > 0
%     set(gca, 'YLim', [0 ylims(2)]);
% else
%     set(gca, 'YLim', [ylims(1) 0]);
% end
 

% ================================== SUBFUNCTIONS =======================================
function handles = analyzeROI(handles)

mapStack = handles.analysis.mapAverager.mapStack;
option = handles.analysis.mapAverager.ROI.option;
switch option
    case 'col'
        rowStart = handles.analysis.mapAverager.ROI.colROIfirstRow;
        rowEnd = handles.analysis.mapAverager.ROI.colROIlastRow;
        colStart = handles.analysis.mapAverager.ROI.colROIfirstCol;
        colEnd = handles.analysis.mapAverager.ROI.colROIlastCol;
    case 'row'
        rowStart = handles.analysis.mapAverager.ROI.rowROIfirstRow;
        rowEnd = handles.analysis.mapAverager.ROI.rowROIlastRow;
        colStart = handles.analysis.mapAverager.ROI.rowROIfirstCol;
        colEnd = handles.analysis.mapAverager.ROI.rowROIlastCol;
    case 'all'
        rowStart = handles.analysis.mapAverager.ROI.allROIfirstRow;
        rowEnd = handles.analysis.mapAverager.ROI.allROIlastRow;
        colStart = handles.analysis.mapAverager.ROI.allROIfirstCol;
        colEnd = handles.analysis.mapAverager.ROI.allROIlastCol;
    otherwise
end
roi = [rowStart rowEnd colStart colEnd];
% if strcmp('col', option)
%     Xvec = handles.analysis.mapAverager.xdata;
% elseif strcmp('row', option)
%     Xvec = handles.analysis.mapAverager.xdata;
% else
%     Xvec = [];
% end   
Xvec = [];


M = mapStack;
[r,c,p] = size(M);

% error catching for mapStack
if p < 2
    disp('mapROImeans: mapStack contains less than 2 maps; cannot evaluate means.');
    return
end
if sum(sum(sum(isinf(M))))>0
    M(isinf(M))=NaN;
end

% error catching for roi
if isempty(roi)
    roi = [1 size(M,1) 1 size(M,2)];
elseif numel(roi) ~=4
    disp('mapROImeans: ROI parameters incorrect.');
    return
end

% calculations
switch option
    case 'row'
%         for n = 1:p
%             m = M(roi(1):roi(2), roi(3):roi(4), n);
%             m = m';
%             V(n,:) = nanmean(m);
%         end
        for n = 1:p
            m = M(roi(1):roi(2), roi(3):roi(4), n);
            m = m';
            if roi(3) < roi(4)
                V(n,:) = nanmean(m);
            elseif roi(3) == roi(4)
                V(n,:) = m;
            end
        end
        vecX = roi(1):roi(2);
        vecY = V;
        vecAvg = nanmean(V);
        vecErr = nanstd(V)./sqrt(sum(~isnan(V)));
        titleTxt = 'average along rows';
        
    case 'col'
%         for n = 1:p
%             m = M(roi(1):roi(2), roi(3):roi(4), n);
%             %         m = m';
%             V(n,:) = nanmean(m);
%         end
        for n = 1:p
            m = M(roi(1):roi(2), roi(3):roi(4), n);
            %         m = m';
            if roi(1) < roi(2)
                V(n,:) = nanmean(m);
            elseif roi(1) == roi(2)
                V(n,:) = m;
            end
        end
        vecX = roi(3):roi(4);
        vecY = V;
        vecAvg = nanmean(V);
        vecErr = nanstd(V)./sqrt(sum(~isnan(V)));
        titleTxt = 'average along columns';
        
    case 'all'
        for n = 1:p
            m = M(roi(1):roi(2), roi(3):roi(4), n);
            V(n) = nanmean(nanmean(m));
        end
        vecX = 1;
        vecY = V;
        vecAvg = nanmean(V);
        vecErr = nanstd(V)./sqrt(sum(~isnan(V)));
        titleTxt = 'average over region';
end

if ~isempty(Xvec)
    if strcmp(option, 'row')
        vecX = -Xvec;
    elseif strcmp(option, 'col')
        vecX = Xvec;
    end
end

% put the output into the handles variable
handles.analysis.mapAverager.ROI.vecX = vecX;
handles.analysis.mapAverager.ROI.vecY = vecY;
handles.analysis.mapAverager.ROI.vecAvg = vecAvg;
handles.analysis.mapAverager.ROI.vecY = vecY;
handles.analysis.mapAverager.ROI.vecErr = vecErr;


